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We examine the stability of self-similar solutions for an accelerating relativistic blast wave which 
is generated by a point explosion in an external medium with a steep radial density profile of a 
power- law index > 4.134. These accelerating solutions apply, for example, to the breakout of a 
gamma-ray burst outflow from the boundary of a massive star, as assumed in the popular collapsar 
model. We show that short wavelength perturbations may grow but only by a modest factor < 10. 
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I. INTRODUCTION 



The self-similar solutions of relativistic blast waves are 
of much interest because of their recent applications to 
the study of Gamma-Ray Bursts (GRBs). A sudden re- 
lease of a large amount of energy within a small volume 
results in a strong explosion that drives a relativistic 
shock into the surrounding medium. At late times the 
blast wave approaches a self-similar phase whereby its 
speed and the distribution of the pressure, density, and 
velocity of the gas behind the shock front do not depend 
on the length and time scales of the initial explosion, but 
only on the explosion energy and the properties of the un- 
shocked external medium. The self-similar solutions de- 
scribing this phase have been first studied by Blandford 
and McKee jj] (hereafter, BMK). We list their central 
results, which are relevant to this paper, in section §11. A. 
Note that in the BMK solution the total energy released 
in the explosion E is the only relevant parameter. 

The BMK solution is only valid for k < 4, where k is 
the power-law index of the radial density profile of the 
external medium, i.e., p\ cx r~ k . When k > 4, the sim- 
ilarity variable defined by Blandford and McKee Q is 
no longer appropriate. Even in the range 3 < k < 4, 
the validity of the BMK solution is not justified, because 
the mass contained behind the shock front diverges if the 
density profile of the shocked fluid is described by the 
BMK solution. 

The self-similar solutions for steep density profiles with 
a power-law index k > 4 were derived recently by Best 
& Sari gp. The derivation of these solutions is similar 
to that in the non-relativistic regime. The self-similar 



solutions to a non-relativistic blast wave were discovered 
independently by Sedov [Q , Von Neumann Q , and Tay- 
lor H . The so-called " Sedov- Von Neumann- Taylor blast 
wave" solutions exist only for k < 5, but Waxman & 
Shvarts |7) showed that in the range 3 < k < 5 these solu- 
tions fail to describe the asymptotic flow because the en- 
ergy diverges; instead they found second-type self-similar 
solutions for 3 < k < 5 as well as for k > 5. The new 
class of non-relativistic, self-similar solutions describe the 
flow in a limited spatial region D(t) < r < R(t), where 
R(t) is the shock radius and D(t) coincides with a C+ 
characteristic so that the flow inside the region r < D(t) 
does not affect the flow in the outer self-similar region. 
The self-similar solution has to cross the sonic line into 
the region where the C+ characteristic can not catch-up 
with the shock front. The solution describes a shock ac- 
celerating with a temporal dependence whose power-law 
index is uniquely determined by requiring that the self- 
similar solution cross the sonic line at a singular point. 
Note that in these second-type self-similar solutions the 
total energy released in the explosion E is not a relevant 
parameter. Although the energy in the self-similar part 
of the flow approaches a constant as time diverges, the 
fraction of the explosion energy E carried by the self- 
similar component depends on the details of the initial 
conditions. Thus, contrary to the BMK case, dimensional 
arguments can not be used to determine the power-law 
index of the temporal dependence. Instead, the singular- 
point determines the temporal power-law index. In the 
ultra-relativistic regime, the second-type self-similar so- 
lutions for k > 4 can similarly be obtained by requiring 
that they cross the sonic line at a singular point. Best 
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& Sari [0 found that these self-similar solutions exist for 
k > 5 — ^/3/4 and describe accelerating shock waves. 
However, the properties of the flow in the self-similar re- 
gion, such as the energy and mass contained in the region, 
were not discussed. 

In this paper, we rederive the self-similar solutions of 
ultra-relativistic blast waves for k > 4 using a differ- 
ent self-similar variable and discuss the properties of the 
flow in the self-similar regime. Our main goal is to study 
the stability of these self-similar solutions. The stabil- 
ity of the Waxman-Shvarts self-similar solutions in the 
non-relativistic regime was studied by Sari, Waxman & 
Shvarts [||. They found that shocks accelerating at a 
rate larger than a critical value and corresponding to so- 
lutions that diverge in finite time, are unstable for small 
and intermediate wavenumbers. Shocks that accelerate 
at a rate smaller than the critical rate are stable for most 
wavenumbers. The acceleration rate can be quantified 
by the measure S = RR/R 2 , where the dots denote time 
derivatives and R(t) is the radius of the shock front. This 
measure provides the fractional change of the velocity 
over a characteristic time scale for evolution (R/R). So- 
lutions that diverge in finite time have 8 > I while others 
have 5 < 1. Thus, when shocks accelerate sufficiently fast 
they become unstable. 

In the following sections we study the stability of the 
self-similar solutions of ultra-relativistic blast waves for 
steep density profiles with a power-law index k > 4. The 
self-similar solutions are described in §11. We list the 
BMK solutions for k < 4 in §11. A and derive the self- 
similar solutions for k > 4 in §11. B. In §111, we discuss 
the properties of the self-similar flow and calculate the 
energy and mass contained in the self-similar regime. In 
§IV and §V, we study the stability of the self-similar so- 
lutions. Finally, we summarize our main results in §VI. 



II. SELF-SIMILAR SOLUTIONS 



A. BMK solutions for k < 4 



dn 1 1 d , , , „ 



(3) 



where n' is the density as measured in the laboratory 
frame, 7 and (3 are the Lorentz factor and velocity of the 
fluid, and 
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is the convective derivative. Throughout this paper we 
set the speed of light c to unity. Assuming that the blast 
wave is ultra-relativistic so that the Lorentz factor of the 
shock front V and the shocked fluid 7 are much larger 
than unity, we only search for solutions accurate to the 
lowest order in j~ 2 and T~ 2 . 

The effective thickness of the blast wave is approxi- 
mately R/T 2 , where R is the radius of the shock front. 
Thus an appropriate choice of similarity variable is 



1 



r 2 > 0. 



(5) 



Next we assume that the external medium has a scale- 
free, power-law density profile p\ cx r~ k . Ignoring ra- 
diative losses, the total energy contained in the shocked 
fluid remains constant and so the Lorentz factor of the 
shock front evolves adiabatically as a power law, 



T 2 oc r 



m > —I. 



(6) 



Keeping only terms up to order 0(T 2 t), the shock ra- 
dius is then given by 



R = t 



I 



1 



(7) 



2(m + I)r 2 

A more convenient similarity variable can be defined as 

X = 1 + 2(m + l)e = [1 + 2(m + l)r 2 ] (l - • (8) 

In terms of x, the pressure, velocity, and density in the 
shocked fluid can be written as 



For pedagogical reasons, we first briefly outline the 
derivation of the self-similar solutions of relativistic blast 
waves for k < 4 by Blandford & McKee. For a complete 
derivation, the reader is referred to the original paper [jjj . 

Assuming an ultra-relativistic equation of state, p — 
(I/3)e, where p and e are the pressure and energy den- 
sity measured in the fluid frame, the equations describing 
a relativistic, spherically-symmetric, perfect fluid can be 
written as, 



dt 



{.PI 4 



2&P 
7 dt' 



(1) 



(2) 



p = gWir 2 /(x), 



:r 2 5 (x), 



»' = 2n 1 r 2 h( X ), 



(9) 



(10) 



(11) 



where X > f, wi and ri\ arc the enthalpy and number 
density of the unshocked external medium. We assume 
that the unskocked external medium is cold, so that w\ 
equals the energy density p±. The jump conditions for a 
strong ultra-relativistic shock are satisfied by the bound- 
ary conditions 



/(l) = fl (l) = fc(l) = l. 



(12) 
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For an adiabatic impulsive blast wave, equations (|])-(||) 
admit a simple analytical solution, first derived by BMK 



7 2 = ^ 2 . 9 (£), 



(21) 



f = X 



-(17-4fc)/(12-3fc) 



.9 = X 



for 



h = x -( 7 - 2fe )/(4- fc ) 7 



m = 3 — k > — 1. 



(13) 
(14) 
(15) 

(16) 



B. Self-similar solutions for > 4 

In searching for self-similar solutions for k > 4, we as- 
sume that the Lorentz factor of the shock front still obeys 
a power law, T 2 oc with m < —1. When m < —1, the 
similarity variable x defined in equation flq) (and used by 
Best & Sari [^) could be negative. For convenience we 
will use £, defined in equation (||), instead as our similar- 
ity variable. If at an initial time to, the shock radius is 
Rq and the Lorentz factor of the shock front is Tq, then 
at a later time t, to 0(Y~ 2 t) the shock radius is given by 



R — Rq + t 



1 



1 



2(m+ l)r 2 
We can rewrite this equation as 



to 



2(m+l)rg 



R = t 



1 



1 



2(m + l)r 2 



a, 



(17) 



(18) 



where a is a constant dictated by the initial conditions. 
This equation for R with m < — 1 differs from equation 
(Q) by a constant a. However, we can choose the initial 
time to such that a is equal to zero. This is appropriate 
because of two reasons. First, the self-similar solutions 
are valid at much later times t » to, thus the effect 
of the special choice of to can be ignored. Second, what 
matters in the derivation of the self-similar solutions is 
the derivative of R, instead of R itself. When a = 0, the 
similarity variable becomes 



l 



i 



2(m+ 1) 



(19) 



Note that we have ignored higher order terms in T~ 2 in 
the above expression. 

Similarly to equations jiil)-(|lll), we write the pressure, 
velocity, and density in the shocked fluid as 



p= 3^r 2 /(0, 



(20) 



n' = 2 ni T 2 h{£,), 
where £ > and the boundary conditions, 
/(0) - g(0) = h(Q) = 1, 



(22) 



(23) 



correspond to the jump conditions for a strong ultra- 
relativistic shock. 

We can now treat T 2 and £ as two new independent 
variables in place of r and t, and get 



d 



d 

dr 2 



1 m 



2 m+ 1 



(m+l)£ 



d_ 



dr 



d 2 d 

t— = —ml -r— t: 
dt dT 2 



2(m + 1) 



d_ 



(24) 



(25) 



~ + (m + !)£-- 
2 g 



w (26) 



In deriving the above equations, we have assumed that 
the blast wave is ultra-relativistic so that T >> 1 and 
7 >> 1. Thus we only keep terms of the lowest contri- 
bution order in T~ 2 and 7~ 2 . 

Substituting equations |24|)-(]26]) into equations (Q)- 
(||), we obtain the following differential equations for /, 
<7, and h: 



2(3m + k)g + [g + 2(m + 1)<?£ + 2]--^ 
+2[<7 + 2(m + l) 9 £-2]-^| =0, 



(27) 



2(5m + 3fc - 8)3 + 3[g + 2(m + 1)#£ - 2] 
+2[ 5 + 2(m+l).g£ + 2]i^| = 0, 



idf_ 



(28) 



2(m + k - 2)g + [g + 2(m + l)g£ - 2]-^ 

1 dg 
+2--| = 0. 

g d£ 



(29) 



Using a new variable, 

y=[l + 2(m + mg, (30) 
we rewrite equations (p^)-(p^) as follows 
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1 

9 



'1 df 



1 / 1 dg 
9 \gd4 



2[4(2m + k-2)-(m + k- 4)y] 

y 2 - 8y + 4 : 



2[(7m + 3fc-4) - (m + 2)y] 



(31) 



(32) 



\ (if)=-2[(2/ 2 -8 2/ + 4)(y-2)]- 1 

x [(m + k - 2)y 2 - (10m + 8fc - 12)y 

+ (18m + lOfc- 16)]. (33) 

One solution to the above equations is obtained for y = 1 
and m = 3 — fc, and corresponds to the BMK solution. 
From the definition of y, y = [l + 2(m+ 1)£]<7, and the re- 
quirement that g be positive it follows that this solution 
is only valid for m > — 1, i.e., k < 4. 

In our search for possible solutions with fc > 4, we start 
by analyzing equations (|3l|)-(|33|). The right-hand-side of 
these equations diverges to infinity if y 2 — 8y+4 = 0. This 
corresponds to two singular points, y\ = 4 — 2-\/3 = 0.536 
and j/2 = 4 + 2V3 = 7.464. In addition, equation j3^ ) 
has another singular point at 2/3 — 2. The solution to 
equation (^) can bypass the singular points j/i and 2/2 
if the numerator on the right-hand-side of the equation 
vanishes at y\ or j/2- This gives 



mi 2 



- 4fc + (fc - 4)2/1,2 
8 - 2/1,2 



(34) 



It is easy to prove that when m = m±, equations ( p2[ ) and 
( |33| ) will also bypass the singular point 2/1 (the numerator 
in the right hand side of each equation is equal to zero 
at 2/1 )■ The same is true for m — m 2 and the singular 
point 2/2- We will show below that when m = mi and 
fc is bigger than a critical value fc c , we have y < 1, thus 
equations d3l"|)-(|33|) are able to bypass the singular point 
2/i and never reach 2/2 and 2/3 • The critical value fc c can 
be calculated by setting mi equal to 3 — fc, the m value 
corresponding to a BMK solution. For y\ = 4 — 2\/3, we 
get 



mi 



12V3- 20+ (3 - 2V3)fc. 



Thus mi = 3 — fc gives us 



k c = 5- — = 4.134 
2 



(35) 



(36) 



The value of mi corresponding to fc c is mi c = —2 + 
v3/2 = —1.134. Thus when fc > fc c , we have mi < mi c . 
We now prove that when m = mi and fc > fc c , we al- 
and the definition 



ways have y < 1. Using equation (|3 
of y in equation (BO) we obtain 



1)5(0 



2/ rfff 



3(0 # 

23(0 [?/ 2 + (m - 3fc + 12)2/ - 4(m + 1)] 

y 2 - 82/ + 4 • 



(37) 



When m = mi, the above equation can be rewritten as 

My - b) 



dy 
d? 



3 - 2/2 



where 



6 = 4- 10V3 + 2\/3fc = 1 + 2\/3(fc - fc c ). 



(38) 



(39) 



When fc > k c , we have 6 > 1. We also have g > and 
2/2 = 4 + 2\/3 > 1, and so the right hand side of equa- 
tion d38| ) is negative when y < 1. Since the boundary 
condition is j/(£ = 0) = 1, y(£) must be a monotonically 
decreasing function of £ with 2/(0 < 1. The asymptotic 
behavior of 2/(0 can be derived as follows. When £_is 
large, y is negative and \y\ is large so that equation ( |3q ) 
can be approximated as 



-23 



1 



mi + U' 



where we have used the approximation y ~ 2 (mi 
for large £. Equation (|o|) yields 

-2/cxr 1/(mi+1) - 



(40) 
1K3 



(41) 



Note that when fc > fc c , we have mi < m\ c < — 1. Thus 
the exponent in the above power law is always positive. 

It can be proven that when fc < k c , equations (|3"l|)- (|33|) 
can not bypass all the singular points with either mi or 
m.2, and so the BMK solution is the only possible solu- 
tion. When fc > fc c , the equations can not bypass all the 
singular points with m2 . But this by itself is not sufficient 
for justifying that mi is the only viable choice. What if 
the solutions cut-off at some radius before reaching any 
singular points? We know that in order not to run into 
divergences of the energy or mass of the system, the solu- 
tions must be truncated at some radius r + (or £ + in terms 
of the similarity variable), which should coincide with a 
C+ characteristic line. The C+ characteristic guarantees 
that the flow in the inner region r < r+ will not influ- 
ence the flow in the self-similar region r + < r < R. This 
C+ characteristic should not overtake the shock front in 
finite time, otherwise the self-similar region will eventu- 
ally disappear. This argument has been applied in the 
non-relativistic case @. We will prove below that in or- 
der to get to the regime where the C+ characteristics can 
not catch the shock front, the solutions have to pass the 
singular point 2/1, making mi the only viable choice. 

First, let us derive the equation for a C + characteristic. 
We use v to denote the fluid velocity in the laboratory 
frame. The sound speed in the fluid frame is u' s = 1/V3- 
Thus the sound speed in the laboratory frame u s is given 
by 

< + « , V3-1 1 
= 1 — 



= 1 - 



„J a ' V3 + l2 7 2 
V3- 1 1 



v's + irw 



(42) 
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where we only keep the first-order term in 7 2 . Thus, a 
C+ characteristic is described by 



dr + 
~~dt 



= 1 - 



V3 - 1 



1 



V3 + ir2 5 (£ + ) 



(43) 



We can rewrite this equation in terms of the similarity 
variable £ + . Using the definition of £ in equation (^J), and 

1 - i/(2r 2 ), dr 2 /dt = -™r 2 /t, 



the relations, dR/dt 
we obtain 

dr _ 1 
dt ~ l ~ 2T 2 



(m 



11c ±^L 

r 2 { r 2 dt ' 



where we only keep the first-order term in T 2 . 
tuting equation into equation (pTJ), we get 



d? H 



\/3-l 1 , 1 



(44) 
Substi- 

(45) 



which describes the evolution of a C+ characteristic. We 
can further rewrite this equation as 



f 



dt 



1 



2g(e+) 



(y+ - yi), 



where 



y + =[l + 2(m+l)S+]g(S+). 



(46) 



(47) 



Equation ( [46] ) implies that when y + > ?/i , the right-hand- 
side of the equation is negative and so £+ will decrease 
with time and the C+ characteristic will approach the 
shock front. Only when y + < y\, £ + will increase with 
time and the C+ characteristic will not overtake the shock 
front. Also notice that the self-similar solution has the 
boundary condition y(£ = 0) = 1 > y\. We thus proved 
that in order to get to the regime where C+ character- 
istics can not overtake the shock front, the self-similar 
solution must pass through the singular point y\, and 
therefore mi is the only viable choice. 

We can now attempt to obtain the self-similar solutions 
for equations (|3l|)-(^3|). For m = mi, these equations be- 
come 



1 df 



1 fldg 
9 \gd£, 



1 /l dh 
~9 \hd£ 



2(mi + k - 4) 
y-V2 



2(mi +2) 

y- 2/2 



2 (mi +k- 2)(y - d) 

(v-V2)(y-2) ' 



where 



d = 4 + V3- 



3 + 2V3 



(48) 



(49) 



(50) 



(51) 



Treating y as the independent variable instead of £ and 
making use of equation (|3S|), equations (^8|)-(50) can be 
rewritten as 



Id/ _ mi + k - 4 
Jdy~~ 



ldg 

gdy 



V 



mi 



y-b 



ldh _ {m 1 + k-2){y-d) 
hdjj ~ (y-2)(y-b) 



(52) 



(53) 



(54) 



The boundary conditions are f(y = 1) = g(y = 1) = 
= 1) = 1. Equations ( j52| ) and ( p3[ ) have the solu- 
tions 



/ = 



3 = 



6- y 



6-1 



mi+fe— 4 



rni+2 



(55) 



(56) 



A special case is obtained at k = 6 (mi = —2) for which 
f(y) =j(y) = l- When 6 ^ 2 (fc ^ (15 - s/3)/3), equa- 
tion (BJ) has the solution 



y 



d-2 / r \ d-b' 

b-y 
6-1 



(mi+/c-2)/(2-b) 



When b = 2 (k = (15 - \/3)/3), equation | 
solution 



(2-W) 



x exp 



mi 2 



(mi + fe- 2)(d- 2) 



- 1 



■ (57) 
has the 

■ (58) 



In general, the functions /(£), g(£) and h(£) do not ad- 
mit simple analytical forms. Their values can be derived 
numerically from equations (p5|)-(p8|). For example, <?(£) 
satisfies the implicit algebraic equation 



9(0 



= / 6-[l + 2(mi + lK]g(Q \ 



mi+2 



1 



(59) 



But generally, we can derive the analytical forms for the 
asymptotic behaviors of /(£), and h(£) in the limit 
of large £. In this limit, y is negative and \y\ is large, and 
so equation (KS) yields 



-.'/ 



mj+2 



-2(mi + 1)^.9 



mi+2 



(60) 



We can solve <?(£) from the above equation and get 

5 (0^a 5 r (mi+2)/(mi+1) , (61) 
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where 



a a 



-2(mi + 1) 
6-1 



-(mi+2)/(mi + l) 



(62) 



Using equation ( |6l| ) we obtain the asymptotic form of 
for large £ 



y(C)^-a,r 1/(mi+1) , 



where 



-2(m! + 1) 



1 -l/(mi + l) 



(6 - l)»U+2 



(63) 



(64) 



Using equation (j63|), we can derive the asymptotic form 
of /(C) f° r large £ from equation (|5q) and get 



/(£) ~ a .£-( m i+ fe - 4 )/( m i+ 1 ) ; 



where 



a/ 



-2(mi + 1) 



6-1 



-(mi+fc-4)/(rai+l) 



(65) 



(66) 



Similarly, the asymptotic form of h(£) for large £ can be 
derived from equations (j57|) and (pq). We obtain 



ft(0=ia h r (roi+fc " a)/(mi+1) . 



(67) 



where 



ah 



-2 (mi + 1) 



1 - (6- l)(mi+l)(d-2)/(2-b) 



-(mi+fe-2)/(m 1 + l) 



(68) 



if 6 ^ 2, and 



-2(m x + 1) 



-(mi+fc-2)/(mi+l) 



(6- l)™i+ 2 _ 

x exp[-(mi + k-2)(d - 2)] 



(69) 



if 6 = 2. 

Equation (^) , which describes the evolution of the C+ 
characteristic, can also be rewritten using the variable y. 
Using equation (38), we obtain 



dy+ = (y+ - vi)(y+ - b) 
dt y+~V2 



(70) 



We have proven earlier that when k > k c one finds 6 > 1 
and y < 1. Thus the sign of the right-hand-side of the 
above equation is decided by the term (y + — y±). If a 
C + characteristic emerges from the region y < yi, i.e., 
y+(t = to) < yi, we know that y+(t) will decrease with 
time or equivalently £ + (t) will increase with time and the 
C+ characteristic will not catch-up with the shock front. 

Equations (|4^)-(|50|) can also be solved numerically for 
different values of k. We plot the results for k = 5.5 



(mi = -1.77) and k = 6.5 (mi = -2.23) in Figure 1. 
When k c < k < 6, the function /(£) decreases with in- 
creasing £ while the function g(£) increases with increas- 
ing £. This implies that when moving inwards away from 
the shock front, the pressure decreases while the Lorentz 
factor increases. For k > 6 the situation is reversed. For 
all k > k c , the function h(£) increases with £ implying 
that the density always increases when moving inwards 
away from the shock front. 



k=5.5, m=-1.77 



k=6.5, m=-2.23 




5 10 5 

FIG. 1. Distributions of the self-similar functions /(£), g(£) 
and 6(f). The left column corresponds to k — 5.5 and 
mi = —1.77, while the right column corresponds to k — 6.5 
and mi = —2.23. 



III. PROPERTIES OF THE FLOW IN THE 
SELF-SIMILAR REGION 



We now examine the properties of the flow in the self- 
similar region bounded by £ = and £ — £+(i), where 
£+(i) coincides with a C+ characteristic that emerges 
from the region y < y±. The evolution of the C+ char- 
acteristic is described by equation (^) with the ini- 
tial condition £+(< = to) = £o, and correspondingly 
y + (t = t ) = y < yi- By solving equation @, we 
get the following equation for y+(t), 



m - y± 
yi - yo 



(,V2-yx) I '(6-1/1) 



(b-y 2 )/(b- Vl ) 



b- yo ) 



t 

to" 

(71) 
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We can now derive the asymptotic behavior of the C+ 
characteristic at t — * oo. In this limit |y+(£)| is large, 
and equation ([7l|) yields 



V+(t) ^ -aot, 



(72) 



where 



= (2/1 - y a ) {y2 - yi)/{b - Vl \b - yo) (b - y2)/{b ' Vl) /h. (73) 



Substituting equation ( p3|) into equation (72), we get 

e+(i) ~ a+f- (mi+1) , (74) 



where 



a + 



-(mi + l) 



The initial condition for the C+ characteristic, namely 
the value of the C+ characteristic with which £+(t) co- 
incides, contains the information about the initial explo- 
sion. 

The shock front is accelerating with a power-law tem- 
poral dependence of its Lorentz factor, T 2 oc t~ mi . How 
does the C+ characteristic propagate? From equations 
( |6l| ) and (|7J) we find that when t — ► 00, 



7 2 (e + ) = ir 2 3 (e + )(xt 2 . 



(76) 



We see that irrespective of the value of k, the C+ char- 
acteristic always accelerates as j 2 oc t 2 . We can also 
calculate the thickness of the self-similar region. Using 
equations (J5J) and (f74j) we obtain when t — ► 00, 



7? 



r 2 



const. 



(77) 



Note that when —2 < mi < m lc , the C + character- 
istic accelerates faster than the shock front, but be- 
cause the Lorentz factors of both surfaces are acceler- 
ating as power-laws of time, the C + characteristic can 
never catch-up with the shock front. Instead, the dis- 
tance between the two surfaces approaches a constant 
value at late times. 

We can now examine the energy and mass contained 
in the self-similar region. The energy contained in the 
spherical shell between £ = and £ = £+(f) is given by 



E 



16npj 2 r 2 dr 



16tt 



Wl T 2 R 3 



(78) 



Using equations ( p5[) and (|6l|), we can calculate the above 
integral for large values of £ + . This gives, 



E 



16tt 



mi + 1 
mi + k — 3 

(mi+k-3)/(mi+l) 



a/a g u>ir 2 i? 3 



(79) 



In deriving the above result we have used the fact that 
when k > k c , we have mi + k — 3 > and mi + 1 < 
0, so that the exponent of in the above equation 
— (mi + k — 3)/(mi + 1) is positive. When t — > 00, 



£+ is given by equation ([74[). 
pi oc Thus when i 



In addition T 2 oc t 



U'l 



E 



const. 



(80) 



The total number of particles contained between £ 
and £ = £+ (t) is given by 



N 



n'4irr 2 dr = 87rni_R 3 



€+ 



(81) 



(75) Using equations (|67j) and J74|), we obtain that for t — > 00, 



fc-3 



const. 



(82) 



We have thus proven that both the energy and mass 
contained between the C+ characteristic and the shock 
front will approach constant values as t — > 00. The situ- 
ation is similar to the non-relativistic case H. 



IV. APPROXIMATE (ANALYTIC) STABILITY 
ANALYSIS 



In order to analyze the stability of the self-similar so- 
lutions obtained in §11. B, we first follow an analytic ap- 
proach (in §IV) based on the assumptions of variable 
separation and a fixed To, where T$ is the unperturbed 
Lorentz factor of the shock front. As we will explain later, 
these assumptions limit the generality of the results. We 
then use numerical simulations (in §V) to directly solve 
the evolution of the perturbations without those assump- 
tions. The numerical simulations demonstrate that the 
results obtained using the analytic approach are qualita- 
tively valid. 



A. Derivation of linear perturbation equations 

For the analytic approach to the stability analysis of 
the self-similar solutions, we use linear perturbation anal- 
ysis similar to that used in the non-relativistic case ||. 
We start from the equations of motion for an ideal rela- 
tivistic fluid: 



d_ 

at 



(n') + V • (n'v) = 0, 



(83) 



7 2 (e +P) 



dv 
~dt 



+ (v • V)v 



Vp + V^)=H. (S4I 
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d_ 

dt 



77/ 4 /3 



+ (V ■ V) 



7 4 / 3 p 

77/ 4 /3 



0. 



(85) 



where v and 7 are the fluid velocity and Lorentz fac- 
tor in the perturbed solution measured in the laboratory 
frame, e and p are the energy density and pressure in the 
perturbed solution measured in the fluid frame, and ri 
is the fluid density in the perturbed solution measured 
in the laboratory frame. We use the Eulerian perturba- 
tion approach, i.e., the perturbed quantities are defined 
as the difference between the perturbed solution and the 
unperturbed one in the same spatial point. Therefore we 
define the perturbed hydrodynamic quantities as 



Sp(r, 9, <p, t) = p(r, 9, <j), t) - p Q (r, t), 



(86) 



27n7 dr 



7o /3 p /2<5 7 2 dp 4, Sri" 



L n 4/3 V3 7 2 Po 3 77^ 



1 (1 1 M 9 K 3po 

"2 7o 4 I + 2 7o V dr { '4/3 



where we have used the relations 



8vq 1 
~dt = H 



dr 



H = 0, 



1 \ H 



2 7 4 ^ + 2^ dr ' 



(92) 



(93) 



(94) 



<5v(r, 9, <j), t) = v(r, 0, 0, £) - 7j (r, t)r = 5w r r + Svt, 



(87) 



<S»'(r,e,0,t)=n / (r,fl,^t)-nJ(r,t), 



where the quantities with subscript "0" are the unper- 
turbed values. Substituting the above quantities into 
equations (|83|)-(|S5|), we obtain the following linear per- 
turbation equations, 



8_ 

at 



Sri 



1 d 



27o 2 



f 

2^ 



1 



f 

27o 



Sri 



I + ?7qVt • Sv T = 0, 



7o 4 



9lo gTo 
2 7 2 7 dt dr 



(lo$P + PoSj 2 ) 



f - 
f 

f 

"2tI 



1 



d 



2 7 2 / <9r 



1 



1 



<5p - 



1 

2tI 

-5p 



1 



SY 



dr J 



(89) 



(90) 



<57J r 



(95) 



and the operator V T = (9/r)(d/d9) + ((f)/r){d/d4>) acts 
as follows on a scalar ^ and a vector f: 



Vt* = - 



f a* 



f a* 



r 96* 



r sin 



Id 1 
V T • f = — r-TT^(sin^) + — 



• sin 6* 96> 



r sin ( 



(96) 



(97) 



Since the unperturbed quantities satisfy equations 
we write 



Pa = gPir§/(0, 



7o = 2 r o5(0, 



n{, = 277 1 r 2 /7(£), 



(98) 



(99) 



(100) 



where To is the unperturbed Lorentz factor of the shock 
front, and £ is the similarity variable defined as 



(101) 



where Ro is the unperturbed radius of the shock front. 
We further define the perturbation variables as 



4po7o 



8 s 

dt SVT 



1 



1 



d 



1 - TT^r -7rSv T 



1 



-<5V7 



2 7 2 J dr 
VrSp - 



dp \ 

-at) Svt 



0, (91) 



<5 7 2 (r, 9, 0, t) = h? 2 5g(OY lm (6, <j>)X{t), (102) 



8v T (r, 9, <P, t) = -JL5v T (OV T Y lm (9, <t>)X{t), (103) 

1 n 



4/3 

7c/ Po 

/4/3 



2^ 

3^T 



<5p 



4<Jn' 

3 77^ 



Sp(r, 9, t) = - Pl r 2 Sf(C)Y lm (9, 4>)X(t), (104) 
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Sn'(r, 8, 0, t) = 2 ni T 2 Sh(OY lm (e, 4>)X(t), (105) 

where the operator Vt = 6(d/d6) + 0(1/ sm6)(d/d<p). 
Note that the variables £ and i are separated in above 
definitions of the perturbations, and so we consider only 
"global" perturbations §,§. The function X (t) measures 
the amplitude of the perturbation relative to the unper- 
turbed values. 

Substituting equations (^)- (|lOC|) and (|l02|)-(|l05|) into 
equations (89)-(B3), we obtain 



mi 



2(mi + 2) 

(y - V2) 



Sh (y - 2) / 1 1 dSh 
~h 2 \g h~dl[ 



-4(mi + 2) 2 (mi + fc - 2)(y - d) 

(y-2/2) (y-y 2 )(y-2) 
11 dSg^ 1(1 + 1) 
9 9 d£ 



S_9 

9 



r 2 
1 



-<5^t = 0, 



(106) 



2(mi + 2)(y-2) 



Sf 
I 



-k - 3mi + g + 

12/ - 2/2) 
(y + 2) (Hdsn 
2 \9fd4) 

2(mi + 2)(y-4) 2(mi + fc - 4) 



2q- 
(W-2) 



(2/ - 2/2) 
1 1 (% 
.9 3 d£ 



(y ~ 2/2) 



9 



= 0, 



(107) 



2 (mi + q+ 1) - 



2(mi + fe-4) 

(2/ -2/2) 



gSv T 



,dSv T Sf 
- (2/ -2 — tt j- = 0, 



(108) 



2(mi + 2)(y-4) 2(mi + fc-4) 



3 9 3 (2/ -2/2) 
8(m!+fc-2)(j/-rf) 
'3 (y-y 2 )(y-2) 



(2/ ~ 2/2) 
% (y- 2) flldfc 
9 9 d£ 



(mi + fe -4)(y - 2) 

(2/ - 2/2) 
(mi + fc - 2) (3/ - d) 



(2/ - 2/2) 



(y- 2) f 11 dSf 

Sh 

T 



2, n . (11 d8h\ „ 



(109) 



where we have used equations (p8|)-(p0|) and the following 
relations 



dt 



—mi- 



r 2 

t ' 



Or 



2(mi + 1) 



(110) 



(111) 



9| _ 1 
~dt ~ 1 



dp 



r »-2KTT)-^ + 1) 



mi 



dt t 



dni n-i 

"dT ~~ ~~ T 



mi 



2(mi + l)rg 



1 - 



mi 



2(mj + l)rg 



2 (mi 



1 



r 2 
1 



(112) 



(113) 



(114) 



(115) 



Note that in deriving equations ( Jl 06| ) ( |l 09| ) we have 
assumed that there is no perturbation in the external 
medium. Moreover, in order to separate variables, X(t) 
has to be a power law in time, X(t) oc t q , where q defines 
the temporal evolution of the perturbation amplitude. If 
the real part of q is positive then the perturbation grows, 
while if the real part of q is negative then the perturba- 
tion decays. 

In equation ( |106| ), the term 1(1 + 1)/Tq is associated 
with causality, namely the fact that a perturbation can 
only propagate at a speed < c/Tq in the transverse direc- 
tion and hence expand across a maximum opening angle 
of ~ 1/Tq. Since To is a function of time, it is not possi- 
ble to achieve a complete separation of variables for this 
equation in contrast with the non-relativistic case. How- 
ever, for any constant value of To we can still calculate 
the power-law index for the growth of the perturbation, 
q. These results are meaningful if we find q > \m\\, so 
that perturbations grow on a time scale shorter than the 
time scale for changes in To. Therefore, the assumptions 
of variable separation and fixed To limit the generality of 
the results. However, even if we find q < \mi\, we should 
still be able to gain an insight into some qualitative prop- 
erties of the perturbation amplitude evolution. 

Equations (106 )— jT09| ) are a complete set of first-order 
differential equations for Sf, Sg, Sh, Svt- After some 
algebraic manipulations, one may write the equations 
for the first order terms dSf/dt;, dSg/dt;, dSh/d^ and 
dSvx/dt; in the following matrix form 



Svt 



A(q,kJ(l + l)/TlO 



( Sf 

Sg 
Sh 
\ Svt 



(116) 



where A is a 4 x 4 matrix. Note that (y 2 — 8y + 4) = 
(2/ — 3/1 ) (2/ — 2/2). Thus, the solutions for the perturba- 
tion variables must pass the same singular point (or the 







sonic line), y\ = 4 — 2V3, as the unperturbed variables. 
Therefore, the value of q can be found by requiring that 
the solutions pass through the singular point y\. This is 
very similar to the non-relativistic case || . 

In o rder to numerically integrate the differential equa- 
tions ( |l!6| ) and derive q we need to specify the bound- 
ary conditions at the shock front when the shock is per- 
turbed. Since the relativistic jump conditions across the 
shock front must be satisfied, we have 



da 

6g = -jt+2(m + q + l) 



(127) 



Another boundary condition results from the require- 
ment that the tangential velocities must be continuous 
across the shock front, yielding 



5vt 



~v t sr. 



(128) 



P : 



(117) Substituting equations (103) and (124) into this equa- 
tion, we get 



ri = 2r 2 ni, 



(118) 



1. 



(129) 



7 2 = -r 2 

7 2 ' 



Equations (|125| )- (|l27j) and (|l29j) are the four boundary 
(119) conditions necessary to solve the perturbation equations. 



where T is the Lorentz factor of the perturbed shock 
front. By linearizing these boundary conditions with re- 
spect to the perturbed quantities, we find 



B. Numerical results 



Sn' ( ^) 5R = 4r$ni±5R- 2fcrgm^, (121) 
or I at Kq 



(122) 



where SR is the deviation of the perturbed shock radius 
R from th e un pe rtur bed shock radius Rq- In deriving 
equations (12C)— (122), we used the relations p\ oc R~ k 
and 



6T 2 



(123) 



where i5r 2 is the deviation of the square of the perturbed 
shock Lorentz factor T 2 from the square of the unper- 
turbed shock Lorentz factor Tq. 
We now define 



6R(9,<p,t) = r)^RoY lm (9,<p)X(t), 

1 n 



(124) 



where 77 is a scale factor that can have an arbitrary value; 
f or c onvenie nce we set r\ = 1. Su bstit uting equations 
(PI), ^)-(|iO0|), (|l0^) , jlQ4| ) and ( fiu5l) into equations 
(120)-()l22|), we find 



5f = ^+2(m + q+l), 



5h= ^ + 2(m + g+l), 
at, 



(125) 



(126) 



Based on the derivations presented in the previous sub- 
section, we may now examine the stability of different 
modes for different values of k. As a particular example, 
we consider the case of A; — 5.5. We derive q for dif- 
ferent values of the mode wavenumber I by integrating 
the differential equations (116) from the shock front to 
its interior and requiring that the solutions pass through 
the singular point. The results are shown in Figure 2. 
The top panel shows the real component of q (Re[<7]) as a 
function of + l)/Lo, where Tq is treated as a scaling 
factor. As mentioned before, Re[q] determines the growth 
rate of the perturbation. In the bottom panel, we plot 
the imaginary component of q (lva[q\), which provides 
the oscillation frequency of the perturbation. Figure 2 
separates the behavior of q into three different regimes: 

• In the regime of small l/T (0 < l/T < 0.87), q is 
a real number and Re[g] is positive, implying that 
the perturbation grows monotonically in time. The 
value of Re [q] increases as I increases. Note that q 
vanishes in the limit of I ~ + 0. This result can be 
derived analytically by comparing two unperturbed 
spherical solutions with different parameters. 

• In the regime of intermediate I /To (0.87 < I/Tq < 
17), (j is a complex number and Re[q] is positive, 
implying that the perturbation grows while oscil- 
lating. As / increases the real part Re[g] decreases 
while the imaginary part Im[q] increases. Note that 
the transition between the real and imaginary solu- 
tions for q occurs at I/Tq ~ 0.87. This result follows 
from causality. When the wavelength of the pertur- 
bation (~ 1//) is smaller than 1/ro, the maximum 
angular separation of two regions that can interact 
with each other, the perturbation can oscillate. 



10 



• Finally in the regime of large I /To (I /To > 17), q 
is a complex number and Re[q] is negative, imply- 
ing that the perturbation decays while oscillating. 
The value of Re[g] decreases (so the absolute value 
of Re[q] increases) as I increases while the value of 
Im[q] increases as I increases. 



soon afterwards it stops oscillating and saturates. Per- 
turbations with small wavenumbers stay in the regime 
of small I /To. The perturbation grows slowly without 
oscillating at the beginning but soon saturates. 

The above results remain qualitatively the same for all 
values of k > 4.134. 




FIG. 2. Perturbation growth rate, q, as a function of 
y/l(l + l)/r . The upper and lower panels show Re[g] and 
Im[g] respectively. 

The actual evolution of a perturbation is shaped by 
the fact that To increases with time as Tq oc t~ mi . If ini- 
tially the wavenumber of the perturbation is sufficiently 
large so that it is in the regime of large I /To, the per- 
turbation will start to decay while oscillating. As time 
progresses, l/T decreases and so both | Re[g] | and lm.[q] 
decrease, the perturbation decays with slower speed and 
oscillates on longer timescales. As soon as the pertur- 
bation enters the regime of intermediate I /To, it starts 
to grow slowly over time and oscillate on even longer 
timescales. The growth rate slowly increases over time, 
but is always limited by the rather small upper bound, 
~Rc[q] < 0.38. Eventually, the perturbation enters the 
regime of small I /To and grows slowly without oscillat- 
ing. As t increases, the growth rate approaches zero, and 
so the perturbation saturates. Therefore, perturbations 
with large wavenumbers (short wavelengthes) grow when 
1 ^5 I /To ^ 10 only by a modest factor. In the case of in- 
termediate wavenumbers, the perturbation goes through 
the two regimes of intermediate and small I /To- There- 
fore it grows slowly with some initial oscillations, but 



V. NUMERICAL SIMULATIONS 

We have verified the above behavior by a direct inte- 
gration of the partial different equations which determine 
the evolution of the perturbation variables, without as- 
suming separability of th e solu tio ns wi th respect to £ and 
t. Instead of equations (102) - (105), we redefined the 
perturbation variables as 



5 1 2 {r,9A,t) = \Tl5g{tt)Y l 



8w T {r,e,4>,t) 



=j*vr(£,t)V T lW 



8p(r, i 



,4>,t) = - Pl Tl8f{t,t)Y lr , 



(130) 



(131) 



(132) 



(133) 



Equations (|106| ) - ( |109| ) were then replaced by four par- 
tial differential equations (PDEs) for the perturbation 
variables 8f(£,t), 8g(£,t) 8h(£,t) and 8vr(£,t). We then 
solved for the evolution of these perturbation variables 
by numerically integrating the PDEs with appropriate 
initial values. In our numerical simulations, the outer 
boundary (£ = 0) is the shock front where the shock 
jump conditions are assumed to b e satisfied. We can 
still define 8R as in equation (124). Then at the outer 
boundary the perturbation variables satisfy 



8n'{r, 9, </), t) = 2 ni T 2 8h^, t)Y lm (0, 0). 



6f(Q, t) = ^-X(t) + 2(mi + T)X(t) + 2t 



dX(t) 
dt ' 



(134) 



8h{0, t) = ^X(t) + 2(mi + l)X(t) + 2t^^-, (135) 
a£ dt 



Sg(0, t) = %X{t) + 2(m x + l)X(i) + 2i^p, (136) 
at; dt 



Sv T {0,t) =X(t). 



(137) 



The inner boundary is chosen to be sufficiently large so 
as to cover the entire similarity region which is bounded 
by a inner C+ characteristic. This way, the values of 
the perturbation variables at the inner boundary can not 
affect the shock front. 
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The numerical simulations gave us the same behavior 
for the perturbations as the previously mentioned analyt- 
ical results for the growth rate q. In Figure 3, we show the 
evolution of X(t) in a numerical simulation with k = 5.5 
and I /T = 75. We plot X(t) in the range t £ [1, 10] in the 
top panel and t £ [1, 100] in the bottom panel. Note that 
X{t) describes the relative displacement of the perturbed 
shock radius from the unperturbed value and we set its 
initial value to be X(t = 1) = 1.0. From Figure 3 we 
see that X oscillates over time over increasingly longer 
timescales and stops oscillating at late times. Its am- 
plitude first decreases, then slowly increases and finally 
saturates. Overall it grows by a factor of ~ 10. These 
results are consistent with the previous discussion on the 
three regimes for the evolution of the perturbations. 



k=5.5, i/r =75 




k=5.5, l/r =75 




FIG. 3. Evolution of X(t) for k = 5.5 and l/F = 75. The 
upper and lower panels show t £ [1, 10] and t £ [1, 100] re- 
spectively. 



VI. CONCLUSIONS 

We have derived the self-similar solutions for an ultra- 
relativistic blast wave in an external medium with a den- 
sity profile pi oc r~ k and k > 4. The solutions exist for 
k larger than a critical value k c = 4.134. They describe 
the flow in the self-similar region bounded by the shock 
front and a C + characteristic. The shock front acceler- 
ates with Lorentz factor T 2 oc t~ mi and mi < —1.134, 
while the C+ characteristic accelerates with Lorentz fac- 
tor 7 2 oc i 2 . The energy and mass contained inside the 



self-similar region approach constant values as time di- 
verges. 

We have found that at large wavenumbers the pertur- 
bations first decay, then grow slowly over time and even- 
tually saturate. The initial decay and the intermediate 
growth are accompanied by temporal oscillations. These 
small wavelength perturbations grow when 1 < I /To < 10 
with an overall factor of ~ 10. At intermediate wavenum- 
bers, the perturbations first grow slowly and then satu- 
rate. The initial growth is also accompanied by tempo- 
ral oscillations. At small wavenumbers the perturbations 
grow monotonically in time but soon saturate. Our re- 
sults also apply to expanding relativistic jets as long as 
the opening angle of the jet is larger than the inverse of 
its Lorentz factor. 

In the collapsar model of gamma-ray bursts, a colli- 
mated relativistic outflow is generated due to the collapse 
of the core of a massive star. The outflow approaches the 
stellar envelope at a modest semi-relativistic speed but 
is expected to accelerate significantly across the sharp 
density gradient at the surface of the star fie]] . Our re- 
sults indicate that in the breakout phase perturbations 
are close to being stable in spherical symmetry. It is still 
possible, however, that the lateral expansion of the jet at 
breakout would be accompanied by instabilities. These 
instabilities may produce variations in the Lorentz factor 
of the jet needed in the internal shock model. They may 
also be responsible for the complex light curves observed 
in most GRBs. Current numerical simulations (l(J lack 
adequate resolution at the stellar surface to follow the 
shock breakout and confirm the instabilities. We leave 
a detailed study of the instabilities associated with the 
lateral expansion of the jet for future work. 
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